do reconstruct_extrad		/*puts together the various sets of moments*/
set more off			
capture ge dd1=tds==1
cap log close
log using model_rwIDrgREIG.log,replace			 



cap drop bob0 bob1 bob2

ge bob0=bobt==0 & bobs==0
ge bob1=bobt==1 & bobs==1
ge bob2=bobt==2 & bobs==2

replace dd=1 if dd==. & yt==ys
replace dd=0 if dd==.


cap drop age_* agef_*
cap drop  agem1
ge agem1=agetA-1
mkspline age_1 4 age_2 9 age_3 15 age_4 19 age_5 26 age_6 = agem1
mkspline agef_1 5 agef_2 10 agef_3 15 agef_4 20 agef_5 = agetA
*26-30, 31-35, 36-40, 41-45, 46-50, 51-60


replace tds = abs(ys-yt) if tds ==.

cap drop startys
ge startys=(ys-(agesA))				

cap drop noverlap
ge noverlap = min(yt,ys) - startys 



capture program drop listpara		/*loads all the model parameters into a list*/
program define listpara
	version 5.0
	global listpar ""
	while "`1'"~=""	{
		global listpar "$listpar `1'"
		mac shift
	}
end
listpara sa0 sa1 sa2 sb0 sb1 sb2 s0i sti s0ti s0r str s0tr  /*
*/ b3_0 b3_0_1 b3_0_2 b3_0_3 b3_0_4 b3_0_5 b3_0_6 /*
*/ b3_1 b3_1_1 b3_1_2 b3_1_3 b3_1_4 b3_1_5 b3_1_6 b3_2 b3_2_1 b3_2_2 b3_2_3 b3_2_4 b3_2_5 /*
*/ b4_0 b4_1 b4_2 b3_01 b3_02 b3_12  b30_0 b30_1 b30_2   /*
*/ pw1 pw2 pw3 pw4 pw5 pw6 pw7 pw8 pw9 pw10 pw11 pw12 pw13 pw14 pw15 pw16 pw17 pw18 pw19 pw20 pw21 pw22 pw23 pw24 pw25 pw26 pw27 pw28 pw29 pw30 pw31 pw32 pw33 pw34 /*
*/ tw1 tw2 tw3 tw4 tw5 tw6 tw7 tw8 tw9 tw10 tw11 tw12 tw13 tw14 tw15 tw16 tw17 tw18 tw19 tw20 tw21 tw22 tw23 tw24 tw25 tw26 tw27 tw28 tw29 tw30 tw31 tw32 tw33 tw34 /*
*/ tco36 tco39 tco42 tco45 tco48 tco51 





capture program drop droppi		/*cleans from estimated parameters left around by previous runs*/
program define droppi
	version 5.0
	while "`1'"~=""	{
		capture drop `1'
		mac shift
	}
end


droppi $listpar resi  


capture drop br1
capture drop br2
capture drop bc1 
capture drop bc2
capture gen br1=yt<=1991				 
capture gen br2=yt>1991

capture gen bc1=ys<=1991
capture gen bc2=ys>1991



capture program drop nlamft15
program define nlamft15
	version 5.0
	if "`1'"=="?" { 
		global S_1 "$listpar"
		nlinit .08 sa0 sa1 sa2 sb0  sb1 sb2  s0i sti s0ti s0r str s0tr  
		nlinit .04 b3_0 b3_0_1 b3_0_2 b3_0_3 b3_0_4 b3_0_5 b3_0_6 b3_1 b3_1_1 b3_1_2 b3_1_3 b3_1_4 b3_1_5 b3_1_6 b3_2 b3_2_1 b3_2_2 b3_2_3 b3_2_4 b3_2_5 
		nlinit .05 b30_0 b30_1 b30_2  
		nlinit .6 b4_0 b4_1 b4_2
		nlinit 0   b3_01 b3_02 b3_12
		nlinit 1 pw1 pw2 pw3 pw4 pw5 pw6 pw7 pw8 pw9 pw10 pw11 pw12 pw13 pw14 pw15 pw16 pw17 pw18 pw19 pw20 pw21 pw22 pw23 pw24 pw25 pw26 pw27 pw28 pw29 pw30 pw31 pw32 pw33 pw34
		nlinit 1 tw1 tw2 tw3 tw4 tw5 tw6 tw7 tw8 tw9 tw10 tw11 tw12 tw13 tw14 tw15 tw16 tw17 tw18 tw19 tw20 tw21 tw22 tw23 tw24 tw25 tw26 tw27 tw28 tw29 tw30 tw31 tw32 tw33 tw34
		nlinit 1 tco36 tco39 tco42 tco45 tco48 tco51    
		exit
	}
	tempvar pv tv  pl tl  b3 beta30 beta31 beta32 b4 b30 sa sb  tco ID IG RE st sat

	
*CORE PERMANENT

	quietly capture gen double `sa'= ($sa0*bob0 + $sa1*bob1 + $sa2*bob2)
	quietly capture gen double `sb'= ($sb0*bob0 + $sb1*bob1 + $sb2*bob2)

	
	quietly capture gen double `ID'= `sa' + `sb'*min(agetA,agesA)  
	quietly capture gen double `IG'= $s0i + $sti*agetA*agesA + $s0ti*(agetA+agesA) 
	quietly capture gen double `RE'= $s0r + $str*agetA*agesA + $s0tr*(agetA+agesA) 
	
	

	
	quietly capture gen double `pv'= `ID' + `IG' if (bobt==0&bobs==0)									

	quietly replace `pv' = `ID' + `IG' + `RE'      if (bobt==1&bobs==1)|(bobt==2&bobs==2)

	quietly replace `pv'= `IG'	if (bobt==0&bobs>0)
	
	quietly replace `pv'= `IG' + `RE' 	if (bobt==1&bobs==2)


	
*CORE TRANSITORY
	quietly  gen double `tco'=($tco36*dco36+$tco39*dco39+$tco42*dco42+ $tco45*dco45+$tco48*dco48+$tco51*dco51+dco54+dco57) 
	
	quietly capture gen double `b30'= ($b30_0*bob0*`tco' + $b30_1*bob1 + $b30_2*bob2)
	
	quietly capture gen double `tv'	= `b30'	  if dd==1&yt==1980&cohot<=1954 & bobt==bobs		/*places the initial condition of the AR in the right place for each cohort*/
	quietly replace 		   `tv' = `b30'   if dd==1&yt==1982&cohot==1957 & bobt==bobs
	quietly replace 		   `tv' = `b30'   if dd==1&yt==1985&cohot==1960 & bobt==bobs
	quietly replace 		   `tv' = `b30'   if dd==1&yt==1988&cohot==1963 & bobt==bobs
	quietly replace 		   `tv' = `b30'	  if dd==1&yt==1991&cohot==1966 & bobt==bobs
	quietly replace 		   `tv' = `b30'	  if dd==1&yt==1994&cohot==1969 & bobt==bobs
	quietly replace 		   `tv' = `b30'	  if dd==1&yt==1997&cohot==1972 & bobt==bobs
	quietly replace 		   `tv' = `b30'	  if dd==1&yt==2000&cohot==1975 & bobt==bobs
	quietly replace 		   `tv' = `b30'	  if dd==1&yt==2003&cohot==1978 & bobt==bobs
	quietly replace 		   `tv' = `b30'	  if dd==1&yt==2006&cohot==1981 & bobt==bobs
	quietly replace 		   `tv' = `b30'	  if dd==1&yt==2009&cohot==1984 & bobt==bobs

	
	quietly capture gen double `beta30'= $b3_0*exp($b3_0_1*age_1 + $b3_0_2*age_2 + $b3_0_3*age_3 + $b3_0_4*age_4 + $b3_0_5*age_5 + $b3_0_6*age_6)
	quietly capture gen double `beta31'= $b3_1*exp($b3_1_1*age_1 + $b3_1_2*age_2 + $b3_1_3*age_3 + $b3_1_4*age_4 + $b3_1_5*age_5 + $b3_1_6*age_6)
	quietly capture gen double `beta32'= $b3_2*exp($b3_2_1*age_1 + $b3_2_2*age_2 + $b3_2_3*age_3 + $b3_2_4*age_4 + $b3_2_5*age_5 )
	
	quietly capture gen double `b3'= (`beta30'*bob0 + `beta31'*bob1 + `beta32'*bob2)
	quietly capture gen double `b4'= ($b4_0*bob0 + $b4_1*bob1 + $b4_2*bob2)
	
	sort bobt bobs cohot cohos tds yt
	quietly replace `tv'=(`b3'+`tv'[_n-1]*`b4'^2) if dd==1&`tv'==. & bobt==bobs 
	sort bobt bobs cohot cohos  yt ys
	quietly replace `tv'=(`tv'[_n-1]*`b4') if `tv'==. & bobt==bobs
	
	
	
	sort bobt bobs yt ys 

	quietly replace `tv'= $b3_01  * (1-	($b4_0 * $b4_1^tds)^(noverlap))/(1-$b4_0 * $b4_1^tds)	 if `tv'==. & bobt==0 & bobs==1 & yt<=ys & noverlap>0
	quietly replace `tv'= $b3_01  * (1-	($b4_1 * $b4_0^tds)^(noverlap))/(1-$b4_1 * $b4_0^tds)   if `tv'==. & bobt==0 & bobs==1 & yt>ys	& noverlap>0

	
	quietly replace `tv'= $b3_02  * (1-	($b4_0 * $b4_2^tds)^(noverlap))/(1-$b4_0 * $b4_2^tds)	 if `tv'==. & bobt==0 & bobs==2 & yt<=ys & noverlap>0
	quietly replace `tv'= $b3_02  * (1-	($b4_2 * $b4_0^tds)^(noverlap))/(1-$b4_2 * $b4_0^tds)   if `tv'==. & bobt==0 & bobs==2 & yt>ys	 & noverlap>0
	
	
	quietly replace `tv'= $b3_12  * (1-	($b4_1 * $b4_2^tds)^(noverlap))/(1-$b4_1 * $b4_2^tds)	 if `tv'==. & bobt==1 & bobs==2 & yt<=ys & noverlap>0
	quietly replace `tv'= $b3_12  * (1-	($b4_2 * $b4_1^tds)^(noverlap))/(1-$b4_2 * $b4_1^tds)   if `tv'==. & bobt==1 & bobs==2 & yt>ys	& noverlap>0

	
	quietly replace `tv'= 0 if `tv'==.
	
	
	

*LOADINGS
	
	quietly gen double `tl'=(dr0+$tw1*dr1+$tw2*dr2+$tw3*dr3+$tw4*dr4+$tw5*dr5+$tw6*dr6+$tw7*dr7+$tw8*dr8+$tw9*dr9+$tw10*dr10+$tw11*dr11)* ///
				(dc0+$tw1*dc1+$tw2*dc2+$tw3*dc3+$tw4*dc4+$tw5*dc5+$tw6*dc6+$tw7*dc7+$tw8*dc8+$tw9*dc9+$tw10*dc10+$tw11*dc11)*br1*bc1+ ///
				(dr0+$tw1*dr1+$tw2*dr2+$tw3*dr3+$tw4*dr4+$tw5*dr5+$tw6*dr6+$tw7*dr7+$tw8*dr8+$tw9*dr9+$tw10*dr10+$tw11*dr11)* ///
				($tw12*dc12+$tw13*dc13+$tw14*dc14+$tw15*dc15+$tw16*dc16+$tw17*dc17+$tw18*dc18+$tw19*dc19+$tw20*dc20+$tw21*dc21+$tw22*dc22+$tw23*dc23+$tw24*dc24+$tw25*dc25+$tw26*dc26+$tw27*dc27+$tw28*dc28+$tw29*dc29+$tw30*dc30+$tw31*dc31+$tw32*dc32+$tw33*dc33+$tw34*dc34)*br1*bc2+ ///
				($tw12*dc12+$tw13*dc13+$tw14*dc14+$tw15*dc15+$tw16*dc16+$tw17*dc17+$tw18*dc18+$tw19*dc19+$tw20*dc20+$tw21*dc21+$tw22*dc22+$tw23*dc23+$tw24*dc24+$tw25*dc25+$tw26*dc26+$tw27*dc27+$tw28*dc28+$tw29*dc29+$tw30*dc30+$tw31*dc31+$tw32*dc32+$tw33*dc33+$tw34*dc34)* ///
				($tw12*dr12+$tw13*dr13+$tw14*dr14+$tw15*dr15+$tw16*dr16+$tw17*dr17+$tw18*dr18+$tw19*dr19+$tw20*dr20+$tw21*dr21+$tw22*dr22+$tw23*dr23+$tw24*dr24+$tw25*dr25+$tw26*dr26+$tw27*dr27+$tw28*dr28+$tw29*dr29+$tw30*dr30+$tw31*dr31+$tw32*dr32+$tw33*dr33+$tw34*dr34)*br2*bc2 + ///
				($tw12*dr12+$tw13*dr13+$tw14*dr14+$tw15*dr15+$tw16*dr16+$tw17*dr17+$tw18*dr18+$tw19*dr19+$tw20*dr20+$tw21*dr21+$tw22*dr22+$tw23*dr23+$tw24*dr24+$tw25*dr25+$tw26*dr26+$tw27*dr27+$tw28*dr28+$tw29*dr29+$tw30*dr30+$tw31*dr31+$tw32*dr32+$tw33*dr33+$tw34*dr34)* ///
				(dc0+$tw1*dc1+$tw2*dc2+$tw3*dc3+$tw4*dc4+$tw5*dc5+$tw6*dc6+$tw7*dc7+$tw8*dc8+$tw9*dc9+$tw10*dc10+$tw11*dc11)*br2*bc1
	
	
	
	quietly gen double `pl'=(dr0+$pw1*dr1+$pw2*dr2+$pw3*dr3+$pw4*dr4+$pw5*dr5+$pw6*dr6+$pw7*dr7+$pw8*dr8+$pw9*dr9+$pw10*dr10+$pw11*dr11)* ///
				(dc0+$pw1*dc1+$pw2*dc2+$pw3*dc3+$pw4*dc4+$pw5*dc5+$pw6*dc6+$pw7*dc7+$pw8*dc8+$pw9*dc9+$pw10*dc10+$pw11*dc11)*br1*bc1+ ///
				(dr0+$pw1*dr1+$pw2*dr2+$pw3*dr3+$pw4*dr4+$pw5*dr5+$pw6*dr6+$pw7*dr7+$pw8*dr8+$pw9*dr9+$pw10*dr10+$pw11*dr11)* ///
				($pw12*dc12+$pw13*dc13+$pw14*dc14+$pw15*dc15+$pw16*dc16+$pw17*dc17+$pw18*dc18+$pw19*dc19+$pw20*dc20+$pw21*dc21+$pw22*dc22+$pw23*dc23+$pw24*dc24+$pw25*dc25+$pw26*dc26+$pw27*dc27+$pw28*dc28+$pw29*dc29+$pw30*dc30+$pw31*dc31+$pw32*dc32+$pw33*dc33+$pw34*dc34)*br1*bc2+ ///
				($pw12*dc12+$pw13*dc13+$pw14*dc14+$pw15*dc15+$pw16*dc16+$pw17*dc17+$pw18*dc18+$pw19*dc19+$pw20*dc20+$pw21*dc21+$pw22*dc22+$pw23*dc23+$pw24*dc24+$pw25*dc25+$pw26*dc26+$pw27*dc27+$pw28*dc28+$pw29*dc29+$pw30*dc30+$pw31*dc31+$pw32*dc32+$pw33*dc33+$pw34*dc34)* ///
				($pw12*dr12+$pw13*dr13+$pw14*dr14+$pw15*dr15+$pw16*dr16+$pw17*dr17+$pw18*dr18+$pw19*dr19+$pw20*dr20+$pw21*dr21+$pw22*dr22+$pw23*dr23+$pw24*dr24+$pw25*dr25+$pw26*dr26+$pw27*dr27+$pw28*dr28+$pw29*dr29+$pw30*dr30+$pw31*dr31+$pw32*dr32+$pw33*dr33+$pw34*dr34)*br2*bc2 + ///
				($pw12*dr12+$pw13*dr13+$pw14*dr14+$pw15*dr15+$pw16*dr16+$pw17*dr17+$pw18*dr18+$pw19*dr19+$pw20*dr20+$pw21*dr21+$pw22*dr22+$pw23*dr23+$pw24*dr24+$pw25*dr25+$pw26*dr26+$pw27*dr27+$pw28*dr28+$pw29*dr29+$pw30*dr30+$pw31*dr31+$pw32*dr32+$pw33*dr33+$pw34*dr34)* ///
				(dc0+$pw1*dc1+$pw2*dc2+$pw3*dc3+$pw4*dc4+$pw5*dc5+$pw6*dc6+$pw7*dc7+$pw8*dc8+$pw9*dc9+$pw10*dc10+$pw11*dc11)*br2*bc1
				
*EVALUATOR
	quietly replace `1'=`pl'*`pv'+`tl'*`tv'

end 
nl amft15 m   dc* dr* , leave  eps(1e-4)






nlpred resi , resid
cap drop total
nlpred total



do diagncobis   /*computes the correct s.e.*/



mat list result



*********************************PREDICTIONS

*tell stata which is the right Varcov
cap program drop foo2
program define foo2, eclass
	tempname b V
	matrix `b' = e(b)
	matrix `V' = e(V)
	matrix define `V' = vp
	ereturn post `b' `V' 
end

foo2

cap drop rho* 
cap drop IGE* 
cap drop share*
cap drop age
ge age=.

ge rhoI=.
ge rhoI_u=.
ge rhoI_l=.

ge rhoS=.
ge rhoS_u=.
ge rhoS_l=.

ge rhoO=.
ge rhoO_u=.
ge rhoO_l=.

ge IGE=.
ge IGE_u=.
ge IGE_l=.

ge shareIG = .
ge shareIG_u = .
ge shareIG_l = .


global  IG  "_b[s0i] + _b[sti]*`a'*`a' + _b[s0ti]*(`a'+`a')"
global  RE  "_b[s0r] + _b[str]*`a'*`a' + _b[s0tr]*(`a'+`a')"
global  IDS " 0.5*(_b[sa1]  + _b[sb1]*`a') + 0.5*(_b[sa2]  + _b[sb2]*`a')"
global  IDF "_b[sa0] + _b[sat0]*(`a'+`a') + _b[sb0]*`a'"


forv a = 0(1)26 {
	local i = `a' +1
	replace age = `a' + 25 in `i'/`i'

	
	nlcom  (_b[s0i] + _b[sti]*`a'*`a' + _b[s0ti]*(`a'+`a')) /   ///
	((0.5*(_b[sa1] + _b[sb1]*`a') + 0.5*(_b[sa2] + _b[sb2]*`a'))+(_b[s0i] + _b[sti]*`a'*`a' + _b[s0ti]*(`a'+`a'))+(_b[s0r] + _b[str]*`a'*`a' + _b[s0tr]*(`a'+`a')))
	sca bi=trace(r(b))
	mat v=r(V)
	cap drop v1 
	cap drop se
	svmat v , n(v)
	ge se=sqrt(v1)
	qui summ se
	sca se=r(mean)
	drop se
	replace rhoI = bi in `i'/`i'
	replace rhoI_u = bi + 1.96*se in `i'/`i'
	replace rhoI_l = bi - 1.96*se in `i'/`i'

	
	nlcom  (_b[s0r] + _b[str]*`a'*`a' + _b[s0tr]*(`a'+`a')) /  ///
	((0.5*(_b[sa1] + _b[sb1]*`a') + 0.5*(_b[sa2] + _b[sb2]*`a'))+(_b[s0i] + _b[sti]*`a'*`a' + _b[s0ti]*(`a'+`a'))+(_b[s0r] + _b[str]*`a'*`a' + _b[s0tr]*(`a'+`a')))
	sca bi=trace(r(b))
	mat v=r(V)
	cap drop v1 
	cap drop se
	svmat v , n(v)
	ge se=sqrt(v1)
	qui summ se
	sca se=r(mean)
	drop se
	replace rhoO = bi in `i'/`i'
	replace rhoO_u = bi + 1.96*se in `i'/`i'
	replace rhoO_l = bi - 1.96*se in `i'/`i'
	
	
	nlcom  ((_b[s0i] + _b[sti]*`a'*`a' + _b[s0ti]*(`a'+`a')) + (_b[s0r] + _b[str]*`a'*`a' + _b[s0tr]*(`a'+`a'))) /   ///
	((0.5*(_b[sa1] + _b[sb1]*`a') + 0.5*(_b[sa2] + _b[sb2]*`a'))+(_b[s0i] + _b[sti]*`a'*`a' + _b[s0ti]*(`a'+`a'))+(_b[s0r] + _b[str]*`a'*`a' + _b[s0tr]*(`a'+`a')))
	sca bi=trace(r(b))
	mat v=r(V)
	cap drop v1 
	cap drop se
	svmat v , n(v)
	ge se=sqrt(v1)
	qui summ se
	sca se=r(mean)
	drop se
	replace rhoS = bi in `i'/`i'
	replace rhoS_u = bi + 1.96*se in `i'/`i'
	replace rhoS_l = bi - 1.96*se in `i'/`i'

	
	nlcom  (_b[s0i] + _b[sti]*`a'*`a' + _b[s0ti]*(`a'+`a')) / ///
	((_b[sa0]  + _b[sb0]*`a')+(_b[s0i] + _b[sti]*`a'*`a' + _b[s0ti]*(`a'+`a')))
	sca bi=trace(r(b))
	mat v=r(V)
	cap drop v1 
	cap drop se
	svmat v , n(v)
	ge se=sqrt(v1)
	qui summ se
	sca se=r(mean)
	drop se
	replace IGE = bi in `i'/`i'
	replace IGE_u = bi + 1.96*se in `i'/`i'
	replace IGE_l = bi - 1.96*se in `i'/`i'
	
	
	nlcom  (_b[s0i] + _b[sti]*`a'*`a' + _b[s0ti]*(`a'+`a')) /   ///
	((_b[s0r] + _b[str]*`a'*`a' + _b[s0tr]*(`a'+`a'))+(_b[s0i] + _b[sti]*`a'*`a' + _b[s0ti]*(`a'+`a')))
	sca bi=trace(r(b))
	mat v=r(V)
	cap drop v1 
	cap drop se
	svmat v , n(v)
	ge se=sqrt(v1)
	qui summ se
	sca se=r(mean)
	drop se
	replace shareIG = bi in `i'/`i'
	replace shareIG_u = bi + 1.96*se in `i'/`i'
	replace shareIG_l = bi - 1.96*se in `i'/`i'
	

	
}

label var age "Age"
label var rhoI "Intergenerational"
label var rhoS "Overall Sibling"
label var rhoO "Residual Sibling"





graph twoway  ///
			rarea rhoI_u rhoI_l age, color(gs14) || ///
			rarea rhoS_u rhoS_l age , color(gs14)  || ///
			rarea rhoO_u rhoO_l age, color(gs14) || ///
			connect rhoS rhoI  rhoO age, symbol(o + t) lcolor(black black black)  mcolor(black black black) ylab(0(0.1) 0.5, angle(0)) ||	///
			, graphregion(color(white)) legend( label(1 "95% C.I.") order(4 5 6 1)) xtick(25(1)51) saving(rhos_rwIDrgREIG.gph, replace)
 

 

save model_rwIDrgREIG,replace 


global rhoS1 "0"
global rhoI1 "0"
global rhoO1 "0"
global IGE1 "0"
global shareIG1 "0"

global rhoS2 "0"
global rhoI2 "0"
global rhoO2 "0"
global IGE2 "0"
global shareIG2 "0"


global c = 0
global d = 0 



forv a = 0(1)13 {
	
	global c = $c+1
	local b = `a' + 14
	
	global rhoI1 "$rhoI1 + (_b[s0i] + _b[sti]*`a'*`a' + _b[s0ti]*(`a'+`a')) / 	((0.5*(_b[sa1] + _b[sb1]*`a') + 0.5*(_b[sa2] + _b[sb2]*`a'))+(_b[s0i] + _b[sti]*`a'*`a' + _b[s0ti]*(`a'+`a'))+(_b[s0r] + _b[str]*`a'*`a' + _b[s0tr]*(`a'+`a')))"
	global rhoO1 "$rhoO1 + (_b[s0r] + _b[str]*`a'*`a' + _b[s0tr]*(`a'+`a')) / 	((0.5*(_b[sa1] + _b[sb1]*`a') + 0.5*(_b[sa2] + _b[sb2]*`a'))+(_b[s0i] + _b[sti]*`a'*`a' + _b[s0ti]*(`a'+`a'))+(_b[s0r] + _b[str]*`a'*`a' + _b[s0tr]*(`a'+`a')))"
	global IGE1  "$IGE1  + (_b[s0i] + _b[sti]*`a'*`a' + _b[s0ti]*(`a'+`a')) / 		((_b[sa0]  + _b[sb0]*`a')+(_b[s0i] + _b[sti]*`a'*`a' + _b[s0ti]*(`a'+`a')))"
	global shareIG1 "$shareIG1 + (_b[s0i] + _b[sti]*`a'*`a' + _b[s0ti]*(`a'+`a')) / ((_b[s0r] + _b[str]*`a'*`a' + _b[s0tr]*(`a'+`a'))+(_b[s0i] + _b[sti]*`a'*`a' + _b[s0ti]*(`a'+`a')))"

	if `b' < 27 {
		global d = $d+1
		global rhoI2 "$rhoI2 + (_b[s0i] + _b[sti]*`b'*`b' + _b[s0ti]*(`b'+`b')) / 	((0.5*(_b[sa1] + _b[sb1]*`b') + 0.5*(_b[sa2] + _b[sb2]*`b'))+(_b[s0i] + _b[sti]*`b'*`b' + _b[s0ti]*(`b'+`b'))+(_b[s0r] + _b[str]*`b'*`b' + _b[s0tr]*(`b'+`b')))"
		global rhoO2 "$rhoO2 + (_b[s0r] + _b[str]*`b'*`b' + _b[s0tr]*(`b'+`b')) / 	((0.5*(_b[sa1] + _b[sb1]*`b') + 0.5*(_b[sa2] + _b[sb2]*`b'))+(_b[s0i] + _b[sti]*`b'*`b' + _b[s0ti]*(`b'+`b'))+(_b[s0r] + _b[str]*`b'*`b' + _b[s0tr]*(`b'+`b')))"
		global IGE2  "$IGE2  + (_b[s0i] + _b[sti]*`b'*`b' + _b[s0ti]*(`b'+`b')) / 		((_b[sa0]  + _b[sb0]*`b')+(_b[s0i] + _b[sti]*`b'*`b' + _b[s0ti]*(`b'+`b')))"
		global shareIG2 "$shareIG2 + (_b[s0i] + _b[sti]*`b'*`b' + _b[s0ti]*(`b'+`b')) / ((_b[s0r] + _b[str]*`b'*`b' + _b[s0tr]*(`b'+`b'))+(_b[s0i] + _b[sti]*`b'*`b' + _b[s0ti]*(`b'+`b')))"
	}
}




nlcom (corr_ab_I : _b[s0ti] / sqrt(_b[s0i]*_b[sti]))  (corr_ab_R : _b[s0tr] / sqrt(_b[s0r]*_b[str])) ///
	(corr_ab_Tot : (_b[s0ti]+  _b[s0tr]) / sqrt((_b[s0i]+_b[s0r])*(_b[sti]+_b[str])))   



nlcom (rhoI1: ($rhoI1)/$c) (rhoO1: ($rhoO1)/$c) (IGE1: ($IGE1)/$c) (shareIG1: ($shareIG1)/$c)  /*
*/     (rhoI2: ($rhoI2)/$d) (rhoO2: ($rhoO2)/$d) (IGE2: ($IGE2)/$d) (shareIG2: ($shareIG2)/$d), post

 	   


cap program drop foo2
program define foo2, eclass
	tempname b V
	matrix `b' = e(b)
	matrix `V' = e(V)
	ereturn post `b' `V' 
end

foo2



nlcom  (rhoI: (_b[rhoI1]*$c/($c+$d) + _b[rhoI2]*$d/($c+$d))) (rhoS: (_b[rhoI1]*$c/($c+$d)+_b[rhoI2]*$d/($c+$d))  +   (_b[rhoO1]*$c/($c+$d)+_b[rhoO2]*$d/($c+$d)) ) ///
	   (rhoO: (_b[rhoO1]*$c/($c+$d)+_b[rhoO2]*$d/($c+$d)))  (IGE: (_b[IGE1]*$c/($c+$d) + _b[IGE2]*$d/($c+$d)) ) ///
	   (shareIG: (_b[shareIG1]*$c/($c+$d)+_b[shareIG2]*$d/($c+$d))) ///
	   (ratioIG: ((_b[rhoI1]*$c/($c+$d) + _b[rhoI2]*$d/($c+$d)))/((_b[rhoI1]*$c/($c+$d)+_b[rhoI2]*$d/($c+$d)) +(_b[rhoO1]*$c/($c+$d)+_b[rhoO2]*$d/($c+$d)) ))



 
log close
 










